//#include "../DWMRICore/Scalar.h"
//
//#include "vtkImageData.h"
//#include "vtkStructuredPoints.h"
//#include "vtkStructuredPointsWriter.h"
//
//char *pAlphaPathname = NULL;
//char *pColorPathname = NULL;
//char *pVTKPathname = NULL;
//
//CScalar *pAlpha = NULL;
//CScalar *pColor = NULL;
//
//void PrintUsage()
//{
//	printf("Usage: CreateVTKResult -alpha <alpha scalar file> -color <color scalar file> -dst <result vtk file>\n");
//}
//
//int ParseArguments(const int argc, char *argv[])
//{
//	for (int i = 1; i < argc; ++i) {
//		if (strcmp(argv[i], "-alpha") == 0) {
//			pAlphaPathname = argv[++i];
//		} else if (strcmp(argv[i], "-color") == 0) {
//			pColorPathname = argv[++i];
//		} else if (strcmp(argv[i], "-dst") == 0) {
//			pVTKPathname = argv[++i];
//		} else {
//			return -1;
//		}
//	}
//
//	if (pAlphaPathname == NULL || pColorPathname == NULL || pVTKPathname == NULL)
//		return -2;
//
//	return 0;
//}
//
//void CreateVTKResult()
//{
//	int w = pAlpha->GetWidth();
//	int h = pAlpha->GetHeight();
//	int d = pAlpha->GetDepth();
//
//	vtkStructuredPoints *pData = vtkStructuredPoints::New();
//	pData->SetDimensions(w, h, d);
//	pData->SetNumberOfScalarComponents(2);
//	pData->SetScalarTypeToDouble();
//	pData->Update();
//	for (int z = 0; z < d; ++z) {
//		for (int y = 0; y < h; ++y) {
//			for (int x = 0; x < w; ++x) {
//				double *pixels = (double*)(pData->GetScalarPointer(x, y, z));
//				pixels[0] = pAlpha->GetScalarFast(x, y, z);
//				pixels[1] = ((double)z / (double)d * 255.0);
//			}
//		}
//	}
//	pData->Update();
//
//	vtkStructuredPointsWriter *pWriter = vtkStructuredPointsWriter::New();
//	pWriter->SetInput(pData);
//	pWriter->SetFileName(pVTKPathname);
//	pWriter->SetFileTypeToBinary();
//	pWriter->Write();
//	
//	pData->Delete();
//	pWriter->Delete();
//}
//
//int main(int argc, char *argv[])
//{
//	if (ParseArguments(argc, argv) != 0) {
//		PrintUsage();
//		return 0;
//	}
//
//	pAlpha = new CScalar();
//	pAlpha->OpenScalarFile(pAlphaPathname);
//	pColor = new CScalar();
//	pColor->OpenScalarFile(pColorPathname);
//
//	CreateVTKResult();
//
//	SafeDelete(pAlpha);
//	SafeDelete(pColor);
//
//	return 0;
//}

#include "../DWMRICore/Scalar.h"

#include <stdio.h>
#include <omp.h>

char *pFiberPathname = NULL;
char *pMASKPathname = NULL;

CScalar *pMASK = NULL;
float radius = 0.0f;

CPoint3F fiber[200];
CVector3F tangent[200];

void PrintUsage()
{
	printf("Usage: CreateVTKResult -fiber <input fiber file> -radius <radius> -mask <result mask file>\n");
}

int ParseArguments(const int argc, char *argv[])
{
	for (int i = 1; i < argc; ++i) {
		if (strcmp(argv[i], "-fiber") == 0) {
			pFiberPathname = argv[++i];
		} else if (strcmp(argv[i], "-radius") == 0) {
			radius = atof(argv[++i]);
		} else if (strcmp(argv[i], "-mask") == 0) {
			pMASKPathname = argv[++i];
		} else {
			return -1;
		}
	}

	if (pFiberPathname == NULL || radius < 0.0001f || pMASKPathname == NULL)
		return -2;

	return 0;
}

void ReadFiber()
{
	FILE *fp = fopen(pFiberPathname, "r");
	for (int i = 0; i < 200; ++i) {
		float x, y, z;
		fscanf(fp, "%f %f %f", &x, &y, &z);
		fiber[i].m_x = x + 24.5f;
		fiber[i].m_y = y + 24.5f;
		fiber[i].m_z = z + 24.5f;
	}
	fclose(fp);

	tangent[0] = fiber[1] - fiber[0];
	tangent[0].Normalize();
	for (int i = 1; i < 200; ++i) {
		tangent[i] = fiber[i] - fiber[i-1];
		tangent[i].Normalize();
	}
}

//int OnFiberBoundary(CPoint3F &pos)
//{
//	float minDis = 100000.0f;
//	float flag = false;
//	for (int i = 0; i < 200; ++i) {
//		CVector3F d = pos - fiber[i];
//		float dis = d.Length();
//		if (dis < minDis) {
//			minDis = dis;
//			d.Normalize();
//			if (fabs(InnerProduct(d, tangent[i])) < 0.2f) {
//				flag = true;
//			} else {
//				flag = false;
//			}
//		}
//	}
//
//	if (flag == true && minDis < (radius + 0.5f) && minDis > (radius - 2.5f)) {
//		return 1;
//	} else {
//		return 0;
//	}
//}
//
//void CreateFiberMASK()
//{
//	int w = pMASK->GetWidth();
//	int h = pMASK->GetHeight();
//	int d = pMASK->GetDepth();
//	
//	omp_set_num_threads(32);
//
//	for (int z = 0; z < d; ++z) {
//#pragma omp parallel for
//		for (int y = 0; y < h; ++y) {
//			for (int x = 0; x < w; ++x) {
//				CPoint3F pos((float)x/8.0f, (float)y/8.0f, (float)z/8.0f);
//				int flag = OnFiberBoundary(pos);
//				pMASK->SetScalar(x, y, z, flag);
//			}
//		}
//	}
//}
//
//int main(int argc, char *argv[])
//{
//	if (ParseArguments(argc, argv) != 0) {
//		PrintUsage();
//		return 0;
//	}
//
//	ReadFiber();
//
//	pMASK = new CScalar();
//	pMASK->CreateScalar(400, 400, 400);
//
//	CreateFiberMASK();
//
//	pMASK->SaveScalarFile(pMASKPathname);
//	SafeDelete(pMASK);
//
//	return 0;
//}

float OnFiberBoundary(CPoint3F &pos)
{
	float minDis = 100000.0f;
	float flag = false;
	for (int i = 0; i < 200; ++i) {
		CVector3F d = pos - fiber[i];
		float dis = d.Length();
		if (dis < minDis) {
			minDis = dis;
			d.Normalize();
			if (fabs(InnerProduct(d, tangent[i])) < 0.1f) {
				flag = true;
			} else {
				flag = false;
			}
		}
	}

	if (flag == true && minDis < (radius + 0.2f) && minDis > (radius - 0.2f)) {
		return 1.0f - ((minDis-radius)*(minDis-radius)*25.0f);
	} else {
		return 0;
	}
}

int main()
{
	char *raidus_pathname = "Z:\\Downloads\\ground-truth-fibers\\ground-truth-fibers\\fibers-radii.txt";
	float radiuses[20];
	FILE *fp = fopen(raidus_pathname, "r");
	for (int i = 0; i < 20; ++i) {
		int temp;
		fscanf(fp, "%d %f", &temp, &(radiuses[i]));
	}
	fclose(fp);

	char *path = "Z:\\Downloads\\ground-truth-fibers\\ground-truth-fibers\\";

	pMASK = new CScalar();
	pMASK->CreateScalar(400, 400, 400);
	
	omp_set_num_threads(16);

	for (int i = 0; i < 20; ++i) {
		char pathname[512];
		sprintf(pathname, "%sfiber-%02d.txt", path, i+1);
		pFiberPathname = pathname;
		ReadFiber();

		radius = radiuses[i];

		int w = pMASK->GetWidth();
		int h = pMASK->GetHeight();
		int d = pMASK->GetDepth();

		for (int z = 0; z < d; ++z) {
			printf("%d\n", z);
#pragma omp parallel for
			for (int y = 0; y < h; ++y) {
				for (int x = 0; x < w; ++x) {
					CPoint3F pos((float)x/8.0f, (float)y/8.0f, (float)z/8.0f);
					float p = OnFiberBoundary(pos);
					float v = pMASK->GetScalarFast(x, y, z);
					if (p > v)
						pMASK->SetScalar(x, y, z, p);
				}
			}
		}
	}

	pMASK->SaveScalarFile("Z:\\Downloads\\ground-truth-fibers\\ground-truth-fibers\\gt.nrrd");

	SafeDelete(pMASK);

	return 0;
}